Autumn Arctic Pacific Sea Ice Dipole as a Source of Predictability for Subsequent Spring Barents Sea Ice Condition

This study uses observational and reanalysis datasets in 1980–2016 to show a close connection between a boreal autumn sea ice dipole in the Arctic Pacific sector and sea ice anomalies in the Barents Sea (BS) during the following spring. The September–October Arctic Pacific sea ice dipole variations are highly correlated with the subsequent April–May BS sea ice variations (r = 0.71). The strong connection between the regional sea ice variabilities across the Arctic uncovers a new source of predictability for spring BS sea ice prediction at 7-month lead time. A cross-validated linear regression prediction model using the Arctic Pacific sea ice dipole with 7-month lead time is demonstrated to have significant prediction skills with 0.54–0.85 anomaly correlation coefficients. The autumn sea ice dipole, manifested as sea ice retreat in the Beaufort and Chukchi Seas and expansion in the East Siberian and Laptev Seas, is primarily forced by preceding atmospheric shortwave anomalies from late spring to early autumn. The spring BS sea ice increases are mostly driven by an ocean-to-sea ice heat flux reduction in preceding months, associated with reduced horizontal ocean heat transport into the BS. The dynamical linkage between the two regional sea ice anomalies is suggested to involve positive stratospheric polar cap anomalies during autumn and winter, with its center slowly moving toward Greenland. The migration of the stratospheric anomalies is followed in midwinter by a negative North Atlantic Oscillation–like pattern in the troposphere, leading to reduced ocean heat transport into the BS and sea ice extent increase.


Introduction
The loss of Arctic sea ice since the late 1970s has been observed by routine satellite missions (Stroeve and Notz 2018). It is one of the most robust features accompanying the anthropogenic warming in the late twentieth century (Meehl et al. 2007;Serreze et al. 2007) and is projected to exacerbate in the future (Overland and Wang 2013;Stroeve et al. 2012). On top of the apparent decreasing trend, the Arctic sea ice exhibits strong natural variability (Kay et al. 2011;Notz and Marotzke 2012;Stroeve et al. 2012). Multiple mechanisms have been investigated to understand the origin of the sea ice variability. Some studies attributed it to atmospheric dynamical and radiative drivers (e.g., Kay et al. 2008;Graversen et al. 2011;Herbaut et al. 2015;Urrego-Blanco et al. 2019), while others emphasized the roles of oceanic and sea ice processes (e.g., Shimada et al. 2006;Yeager et al. 2015;Zhang 2015;Årthun et al. 2017;Oldenburg et al. 2018) and their interactions (e.g., Nakanowatari et al. 2014;Krikken and Hazeleger 2015). The interplay of natural variability and decreasing trend leads to spatial heterogeneity of the sea ice variability Zhang 2015;Lee et al. 2017;Onarheim and Årthun 2017).
One of the sea ice loss hotspots lies within the Barents Sea (BS). Its sea ice loss contributes the largest portion to the pan-Arctic sea ice declining trend during boreal winter and spring (Onarheim et al. 2015(Onarheim et al. , 2018, and it has experienced strong natural fluctuations at interannual, decadal, and longer time scales in the past decades (Onarheim and Årthun 2017). Many studies showed that the BS sea ice variability has likely exerted profound impacts on local weather and climate (Kohnemann et al. 2017;Pedersen and Christensen 2019), ecosystem (Dalpadado et al. 2014), Arctic commercial shipping (Eicken 2013;Smith and Stephenson 2013;Melia et al. 2016;Pizzolato et al. 2016;Laliberté et al. 2016), fishery activities (Hollowed et al. 2013;Haug et al. 2017), and natural resource extraction (Jung et al. 2016). Extensive studies also argued that BS sea ice change could have far-reaching influences on weather and climate in lower latitudes via altering large-scale atmospheric circulations in a way that resembles the North Atlantic Oscillation (NAO) (e.g., Cohen et al. 2014;Kim et al. 2014;Overland et al. 2015;García-Serrano et al. 2015;Kretschmer et al. 2016;Sorokina et al. 2016;Screen 2017;Zhang et al. 2019; and many others), although this topic has been controversial (e.g., Overland et al. 2016;Screen et al. 2018;Blackport et al. 2019;Cohen et al. 2020;Peings 2019;Blackport and Screen 2020;Dai and Song 2020;Liang et al. 2020). Therefore, improving our understanding of the processes affecting the sea ice variability in the BS and developing an accurate BS sea ice prediction system may benefit environments and socioeconomics both locally and remotely.
Using an operational sea ice prediction system, Bushuk et al. (2017) found that detrended Barents-Kara sea ice prediction can be skillful at 5-11-month lead time in winter, while summer and autumn counterparts is only skillful at lead times less than Denotes content that is immediately available upon publication as open access. 3 months. The local sea ice reemergence associated with preceding subsurface ocean temperature anomalies was identified as a source of winter Barents-Kara sea ice predictability, as reported in other studies (Bushuk et al. 2015;Bushuk and Giannakis 2017;Bushuk et al. 2019a). Using perfect model sea ice prediction systems, some studies reported that the skillful Barents-Kara sea ice prediction can be traced further back to 1.5-2.5 years (e.g., Day et al. 2014;Bushuk et al. 2019b). The difference between operational and perfect model sea ice prediction indicates that there is room for prediction skill improvement. Furthermore, Bushuk et al. (2015) showed that the sea ice anomalies in the Barents-Kara Seas are correlated with those in the Bering Sea when the latter leads the former by a few months to about 1 year. Hence, inclusion of preceding sea ice information in the Arctic Pacific sector may improve the regional sea ice seasonal prediction in the BS.
In this study, we aim to explore a new source in the Arctic Pacific sector of springtime BS sea ice predictability and examine the underlying mechanism using observational and reanalysis datasets during 1980-2016. Data and analysis methods used in this study are described in section 2. In section 3, we identify the preceding sea ice variations in the Arctic Pacific sector, manifested as a dipolar structure, relating to the BS sea ice variability in the following spring, and use an empirical model to address this potential source of prediction skill of the springtime BS sea ice condition. A series of analyses are then performed to determine the drivers of the autumn sea ice dipole and of the spring BS sea ice anomalies. Finally, we investigate the possible linking processes in establishing this regional sea ice relationship. A summary and discussion are provided in section 4.

a. Observational and reanalysis datasets
Monthly mean sea ice concentration (SIC) and sea surface temperature (SST) during 1980-2016 are obtained from the Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) dataset with a horizontal resolution of 18 3 18 (Rayner et al. 2003). We also use the SIC dataset from the National Snow and Ice Data Center (NSIDC; Peng et al. 2013), which provides very similar results as the correlation between corresponding sea ice indices is always higher than 0.95. Hence, we only show results using the HadISST SIC data.
Monthly and daily zonal and meridional winds, geopotential height, total cloud cover, and atmospheric surface heat fluxes, including latent heat, sensible heat, shortwave, and longwave, are obtained from the European Centre for Medium-Range Weather Forecasts (ERA-Interim) from 1980 to 2016 (Dee et al. 2011), regridded to the horizontal resolution of 18 3 18.
The sea ice thickness (SIT), sea ice velocity (zonal and meridional components), sea ice (thickness) advection [= Á (uh)], ocean-to-sea ice heat flux used to melt sea ice, and subsurface oceanic temperature and velocity are obtained from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS; Zhang and Rothrock 2003). Because the temporal coverage of subsurface oceanic fields (from the surface to 60-m depth) ends in 2015, we only analyze the period of 1980-2015 when using PIOMAS reanalysis data. The sea ice advection and ocean-to-sea ice heat flux are obtained in units of meters per second (m s 21 ), representing how much sea ice thickness the heat flux or mass advection can grow or melt (J. Zhang 2020, personal communication). To convert these values to corresponding values in watts per square meter (W m 22 ) in order to more directly compare with air-sea heat fluxes, we multiply the above unit by sea ice density (920 kg m 23 ) and latent heat of sea ice (0.334 J kg 21 ). As the PIOMAS is based on a hindcast simulation of an ocean-sea ice coupled model with satellite-based sea ice concentration data assimilated, the ocean-to-sea ice heat flux is actually a representation of the net surface heat flux out of the ocean at each grid. Strictly speaking, it corresponds to the ocean-to-sea ice heat flux only for the portion of each grid covered by the sea ice plus the ocean-to-atmosphere heat flux for the sea ice-free surface area. Therefore, we estimate the ocean-to-sea ice heat flux by multiplying the SIC by the PIOMAS surface heat flux at each grid point, and the ocean-to-atmosphere heat flux by multiplying it by one minus SIC. For consistency with the surface heat fluxes provided by ERA-Interim, which is positive downward, we change their sign to obtain the sea ice-to-ocean and atmosphere-to-ocean surface heat fluxes from PIOMAS data. It should be noted that the sea ice temperature can be lower than the freezing point, and the amount of sea ice may not immediately decrease or increase when gaining or losing heat fluxes from atmosphere or ocean due to small but finite heat capacity of sea ice. This effect only affects our interpretation minimally. However, the PIOMAS dataset does not provide sea ice temperature or the amount of heat used for phase change of seawater near the freezing point, so accurate estimates cannot be made here.
The monthly climatology in the 1980-2010 period is removed from all fields to retrieve the anomalies. We also remove quadratic trends to isolate interannual variability from the effect of the global warming and nonlinear Arctic sea ice decline . Similar results can be derived based on the variables after removing linear trends. In the following analyses, bimonthly anomalies are considered.

b. Indices
This study uses two sea ice indices to characterize the temporal evolution and strength of the regional SIC anomalies. The Arctic Pacific dipole sea ice index is defined as the areaweighted SIC anomalies averaged over the East Siberian-Laptev Seas (1458-1788E, 688-788N; western magenta box in Fig. 1a) and divided by their standard deviation minus those over the Beaufort-Chukchi Seas (1778-1448W, 688-788N; eastern magenta box in Fig. 1a) similarly divided by their standard deviation. It is noted that the value of sea ice dipole index does not change significantly without standardization. The BS sea ice index is defined as the area-weighted SIC anomalies averaged over the Barents Sea (408-608E, 688-788N; green box in Fig. 1a).
The monthly NAO index is obtained from National Oceanic and Atmospheric Administration's (NOAA's) Climate Prediction Center, and determined as the leading rotated empirical orthogonal function (EOF) mode of geopotential height anomalies at 500-hPa level in the Northern Hemisphere (208-908N) (Barnston and Livezey 1987). All indices are standardized by removing their mean and dividing by their standard deviation before regression analysis.

c. Linear prediction model and cross-validation
To examine the predictability of spring BS sea ice condition, an empirical forecast model is constructed using univariate or multivariate linear regression (Dunstone et al. 2016;Wang et al. 2017). We train the models with the take-N-year-out cross-validation method (Wang et al. 2017) to exclude the predictand and predictor time series from the prediction period. When N 5 12, for example, the predictand and predictor time series are divided into segments with consecutive 12-yr data. We remove the first segment and train the models with the remaining portion of the time series. The resultant regression coefficients are then used to predict the value in the first 12 years. We repeat the same procedure but exclude the successive segments to obtain the cross-validated predicted time series in the whole analysis period. We only show predicted results with N 5 12, but N 5 3 and 6 yield similar results and slightly higher forecast skills. The forecast skill is measured by the anomaly correlation coefficient between the original and the predicted time series.

d. Statistical significance
The statistical significance for the correlations and regressions is assessed based on the two-sided Student's t test. The effect of serial correlation is taken into account by using an effective sample size (ESS) given as where N is the length of time series, and R x and R y are the lag-1 autocorrelations of time series x and y, respectively (Bretherton et al. 1999). To make the estimate more stable and consistent, we calculate the lag-1 autocorrelations using the full data length. The statistical significance of regression maps could be overinterpreted simply based on local statistical tests at each grid point (Livezey and Chen 1983;Wilks 2016). To address this issue, we adopt the ''field significance'' approach (Wilks 2016) that controls the false discovery rate (FDR). We choose a FDR 5 0.1 to achieve global test level a Global 5 0.05, assuming a spatial decorrelation scale of ;1.5 3 10 3 km in the longitudelatitude domain (see

e. Horizontal ocean heat flux
When investigating the physical processes leading to the BS SIC anomalies, we calculate the horizontal ocean heat flux into the BS domain following Lee et al. (2004). Their reformulation of horizontal heat flux aims to separate the effect of internal heat redistribution, which does not contribute to the total heat content, from that of external heat sources or sinks. The zonal heat flux at the western interface of a rectangular box domain at a specific time step t, for example, is written as where C p is the specific heat of seawater (3850 J kg 21 K 21 ), r o is the density of seawater (1025 kg m 23 ), u is the zonal velocity, T is the ocean temperature, T m is the volume-averaged (over the BS domain) ocean temperature, T* 5 T 2 T m , and the hi operator indicates C p r 0 ÐÐ west ( ) dydz. To examine the relative contributions of anomalous circulation and temperature, we decompose the zonal heat flux into four terms: where X means the time-averaged value of variable X, and X 0 5 X 2 X. The four terms correspond to mean temperature difference carried out by mean zonal flow hu T*i, anomalous temperature difference carried out by mean zonal flow huT* 0 i, mean temperature difference carried out by anomalous zonal flow hu 0 T*i, and anomalous temperature difference carried out by anomalous zonal flow hu 0 T* 0 i. Similar formulation can be applied to the meridional heat flux (replacing zonal velocity u by meridional velocity y) at the northern interface of the domain.

f. Removal of tropical Pacific variability
Because the tropical Pacific variability, including El Niño-Southern Oscillation, is one of the dominant sources of variability in the extratropical atmospheric circulations (e.g., Yeh et al. 2018;Domeisen et al. 2019), we remove the teleconnection effects from all fields. The tropical Pacific variabilities are represented by the three leading EOFs of the monthly SST anomalies in the tropical domain (1008E-1008W, 158S-158N), and their influence is removed by linear regression onto the corresponding principal components, assuming that the tropical Pacific modes lead the extratropical atmospheric and sea ice fields by 1 month. The regressions are calculated separately for positive and negative values of the principal components to account for the asymmetric responses to El Niño-type warming and La Niña-type cooling. Note that similar results are obtained without removing the tropical teleconnections. For example, the correlation coefficients between the two sea ice indices are 0.714 and 0.705, respectively, with and without tropical Pacific variability modes removed. This indicates that the tropical Pacific variability does not influence much the relationships discussed in this study.

a. The predictability of the spring Barents Sea sea ice condition
The seasonal evolution of the BS SIC shows that large sea ice extents are formed during boreal late winter and early spring (i.e., January-May), whereas its largest interannual variability (6.5%), measured by the standard deviation, occurs in April-May (AM) (not shown). Hence, we investigate the sea ice variability throughout the Arctic that may lend potential predictability to the AM BS sea ice. The largest remote linkage is found with the Arctic Pacific sea ice in the preceding September-October [hereafter SO (21)]. Lag regression of the bimonthly-averaged SO(21) SIC anomalies over the Arctic domain onto the AM BS sea ice index shows a significant dipolar SIC anomalies in the Arctic Pacific sector (Fig. 1a). The dipolar structure is characterized by positive anomalies in the East Siberian-Laptev Seas and negative anomalies in the Chukchi-Beaufort Seas, indicating that a sea ice increase in the former region and decrease in the later region precede AM BS sea ice increase by as much as 7 months. We also present the SO(21) SIC anomalies regressed onto the SO(21) sea ice dipole index to illustrate the dipolar structure (Fig. 1b). Reversely, AM SIC anomalies regressed onto SO(21) dipole sea ice index show large positive SIC anomalies in the BS (Fig. 1c), confirming that springtime BS SIC anomalies indeed relate to preceding autumn sea ice dipole anomalies. We also find very similar temporal evolutions of the two sea ice indices during 1980-2016 ( Fig. 1d) and their correlation is as high as 0.71, indicating that ;50% of BS sea ice variability can be explained by the SO(21) sea ice dipole variability.
The correlation between the two sea ice indices at various lead-lag times is presented as a correlation matrix in Fig. 2. Indeed, the 7-month connection between BS and Arctic Pacific sea ice dipole variabilities is the strongest (r 5 0.71) as shown by the 7-month lag correlation (cyan filled circle in Fig. 2). The correlations remain significant 2-3 months before and after AM, implying that the development and decay of BS sea ice also relate to the SO(21) Arctic Pacific sea ice dipole. The autocorrelation matrix of AM BS sea ice index (Fig. 2b) shows significant correlations up to 4-5 months before (i.e., the negative lag on the y axis), corresponding to preceding winter months, without overlap with the preceding autumn months when the cross-correlation with the Arctic Pacific sea ice dipole peaks. It also shows that the autumn BS SIC anomalies have limited persistence that does not extend to the following winter and spring months. The persistence of the BS SIC anomalies exhibits strong seasonality and the largest persistence is found starting from December-January for the following 5-6 months.
The high correlation between the AM BS sea ice index and the SO(21) Arctic Pacific dipole sea ice index suggests that part of the predictability for the springtime BS sea ice condition may come from the previous autumn Arctic Pacific sea ice dipole. An empirical forecast model is constructed based on linear regression (section 2c). The prediction of the AM BS sea ice condition based on the SO(21) Arctic Pacific dipole sea ice index alone captures the evolution of the AM BS sea ice index (green line in Fig. 3a) with the overall anomaly correlation coefficient (ACC) between predicted and original BS sea ice index as high as 0.70 during 1980-2016.
To examine how the prediction skill evolves during our analysis period, we calculate the ACC between the predicted BS sea ice index and the original index with a sliding 12-yr window (magenta line in Fig. 3b). Adding the January-February (JF) NAO as a predictor (for reasons discussed in section 3d) increases the ACC, especially during the last 12-yr period (cyan line in Fig. 3b). The ACCs range between 0.54 and 0.85, with the highest values occurring during the first and last few 12-yr periods and the lowest ones from the mid-1990s to mid-2000s. Considering the cross validation, the ACC for the last 12-yr period, 0.84, can be considered as a true prediction skill in a strict sense, since it is based only on the data points from preceding years. Note that using JF NAO as the only predictor strongly reduces the prediction skill (the ACC is about 0.3 over the whole period, and about 0.6 in the first few 12-yr period and much lower afterward).
To inform the spatial distribution of the prediction skill, the ACC is given at each grid point within and nearby the BS (Figs. 3c,d). The high ACCs are located close to the Severny-Yuzhny Islands and Russian coastal regions, becoming weaker westward and northward. The highest ACC within the BS domain is 0.76 (0.77 with JF NAO information included), while only two grid points present negative ACCs, and both are small and not statistically significant at the 5% level. The average of ACC values within the BS domain is 0.40 (0.39 with JF NAO information included), which is lower than the ACC of the BS sea ice index because of uncorrelated variability in the western portion of the BS.

b. The physical interpretation of the autumn Pacific Arctic sea ice dipole anomalies
Although the autumn sea ice dipole anomalies lend predictability for the subsequent spring BS sea ice anomalies, we do not find a very high correlation (r 5 0.46) between the two poles of the SO(21) sea ice dipole. This raises a question whether the sea ice dipole anomalies represent a physical variability pattern or the dipolar structure is simply a statistical artifact, although a similar autumn sea ice dipole was reported in a seasonal sea ice prediction study of Bushuk et al. (2015), and manifested as a SIC footprint in one sea ice reemergence mechanism and likely relating to Arctic atmospheric circulation (see their Fig. 15). In this section, we intend to illustrate that the autumn sea ice dipole anomalies have both statistical and physical meaning, and result from several key physical processes.
To verify whether the SO(21) sea ice dipole anomalies can be identified using standard statistical tools, we perform an EOF analysis of the SO(21) SIC anomalies in the Arctic Pacific sector. The third leading EOF mode (hereafter EOF3), explaining 12% of the total variance, presents the SIC dipolar The negative (positive) number on the y axis represents how many months the Arctic Pacific dipole sea ice index leads (lags) the BS sea ice index. Open circles and filled circles denote that the correlation coefficients are significant at the 5% and 1% levels, respectively. The highest correlation is indicated by cyan filled circle. (b) As in (a), but for lead-lag autocorrelation coefficients for the BS sea ice index. structure and its principal component (hereafter PC3) is highly correlated with the SO(21) SIC dipole index (r 5 0.84) and the AM BS SIC index (r 5 0.70) (not shown). Therefore, the EOF3 can also represent the autumn SIC dipole and the PC3 characterizes its temporal variability. However, the EOF3 is not the leading EOF mode and it is not statistically separable from the second EOF according to North's rule (North et al. 1982). Hence, we prefer to use the SO(21) SIC dipole index. To also demonstrate that dipolar SIC anomalies indeed occur in our analysis period of 1980-2016, we conduct case studies based on the strength of the sea ice dipole index (magenta line Fig. 1c) and present the SO(21) SIC anomalies in 1998 and 1991 (Figs. 4a,c), which represent the sea ice dipolar anomalies in its positive and its negative phase, respectively. Importantly, the two cases are followed by BS sea ice anomalies in AM (Figs. 4b,d), reflecting the strong relationship between the autumn sea ice dipole anomalies and subsequent spring sea ice anomalies in the BS identified in section 3a. Similar positive and negative sea ice dipole cases with various strength are identified in other years. This indicates that the sea ice dipole anomalies not only occasionally emerge during past decades, but also with different phases.
Next, we perform regression analysis onto the SO(21) sea ice dipole index with up to four months lead time [i.e., back to MJ(21)] to examine the evolutions of the physical properties and potential drivers of the sea ice dipole anomalies. Regressed SIC anomalies from MJ(21) to SO(21) (first column in Fig. 5) indicate that the sea ice dipole anomalies start to emerge as early as JJ (21), albeit with weaker magnitude and a shifted spatial pattern. The negative anomalies in the Beaufort-Chukchi Seas (i.e., eastern magenta box) first appear in the Beaufort Sea in MJ(21) and expand westward into the Chukchi Sea, becoming stronger as time evolves, while the positive anomalies emerge in the East Siberian-Laptev Seas (i.e., western magenta box) in JJ(21) and intensify until SO(21). Regressed SIT and SST anomalies (second and third columns in Fig. 5) correspond well with the SIC anomalies evolution. The consistent evolution of SIC, SIT, and SST suggests that the spring-to-autumn SST and SIT sea ice reemergence mechanisms (Bushuk et al. 2015; is not involved in the development of the autumn sea ice dipole. We proceed to analyze the regressed heat fluxes from the atmosphere into the ocean or sea ice (from ERA-Interim) and those from the sea ice into the ocean (from PIOMAS), and sea ice mass convergence by lateral advection (also from PIOMAS), all of which can physically contribute to sea ice growth or retreat. Their regression maps clearly show that the atmospheric heat fluxes (fourth column in Fig. 5) are the dominant driver of the sea ice dipole anomalies. Starting from MJ(21) to JA(21), the negative atmospheric heat flux anomalies over the East Siberian-Laptev Seas indicate less heating from atmosphere, leading to sea ice increase, while the positive heat flux anomalies in the east of the Beaufort Sea melt more sea ice. These spatial patterns correspond well to the SIC dipole anomalies from MJ(21) to JA(21). When the sea ice dipole anomalies become mature, in AS(21), the dipolar atmospheric heat flux anomalies disappear and then flip sign in SO (21), becoming a response to the sea ice dipole anomalies. It is noted that no strong regressed sea ice-to-ocean heat fluxes (fifth column in Fig. 5) occur from MJ (21) to JA(21), indicating almost no heat exchange between sea ice and the ocean. In other words, most of the atmospheric heat fluxes are used to melt or grow the sea ice, rather than penetrating into the ocean. The regressions of the sea ice velocity and advection do not show coherent anomalies that would drive the sea ice dipole (sixth column in Fig. 5), although there are some localized impacts, e.g., the negative advection anomalies (i.e., divergence of sea ice mass) in part of the eastern box from JJ(21) to AS (21).
The decomposition of atmospheric heat fluxes into their turbulent (sensible plus latent), longwave, and shortwave radiative components (Fig. 6) shows that persistent shortwave dipolar anomalies from MJ(21) to SO(21) plays the dominant role in forcing the sea ice dipole anomalies (fourth column in Fig. 6). Strong shortwave dipolar anomalies appear in MJ(21) and JJ(21), and weaken afterward. They are largely related to changes in total cloud cover (fifth column in Fig. 6). More cloud cover over the East Siberian-Laptev Seas reduces incoming solar radiation, and vice versa over the Beaufort Sea in JJ(21). Additionally, the positive sea ice-albedo feedback mechanism (Winton 2006;Serreze and Barry 2011;Screen et al. 2012) seems at play, in which more incoming shortwave melts more sea ice and produces more open, darker waters, resulting in more shortwave absorption in the ocean and more sea ice melting in subsequent months, and vice versa for less incoming shortwave. On the other hand, significant turbulent and longwave heat fluxes do not appear in the sea ice dipole region until JA(21), but they gradually intensify after JA(21) when the sea ice dipole becomes mature. This can be interpreted as a response to the sea ice dipole because less sea ice cover implies more open water that favors heat release into the atmosphere via turbulent and longwave fluxes, and vice versa for more sea ice cover.
In summary, we find that the Arctic Pacific sea ice dipole anomalies indeed occurred throughout the past decades, seemingly primarily driven by anomalous shortwave radiation related to Arctic cloud cover variability.

c. The physical interpretation of the spring Barents Sea sea ice anomalies
We now perform similar regression analyses onto the AM BS sea ice index to investigate the potential drivers of the spring BS sea ice anomalies. Consistent evolutions of SIC, SIT, and SST anomalies (first three columns in Fig. 7) again suggest that SST or sea ice reemergence is not involved. Instead, the AM BS sea ice anomalies are primarily driven by the heat fluxes (fourth and fifth columns in Fig. 7). Increased heat flux since D(21)J from the sea ice into the ocean (positive values in fifth column in Fig. 7), in other words reduced heat used to melt the sea ice, favors BS sea ice growth. On the other hand, there is increased atmospheric heat flux into the BS (positive values in fourth column in Fig. 7), part of which could be used to melt sea ice. The sea ice advection (sixth column in Fig. 7), in contrast, shows a west-east dipolar anomalies within the BS, which neither corresponds well to the overall positive BS sea ice anomalies nor contributes to BS sea ice increase in terms of the BS-averaged anomalies. This implies that the growing anomalous BS sea ice extent during winter to spring likely results from the competition between atmospheric and oceanic heat fluxes used to melt or grow the sea ice.
To investigate the oceanic processes that contribute to the increase of BS sea ice, we next perform a heat budget analysis over the BS ocean domain bounded by land to the south and east, and by 60 m at the bottom, the deepest level that the FIG. 5. Regression maps onto the SO(21) dipole sea ice index from MJ(21) to SO(21) for anomalous (first column) SIC, (second column) SIT, (third column) SST, (fourth column) atmospheric net surface heat flux (positive downward, from ERA-Interim reanalysis dataset), (fifth column) sea ice-to-ocean heat flux (positive downward, from PIOMAS), and (sixth column) sea ice velocity (arrows, from PIOMAS dataset) and convergence of sea ice via lateral advection (from PIOMAS dataset). The regression coefficients within the magenta contour lines are significant at the 5% level, while the cyan contour lines inform the field significance. The regression coefficients of sea ice velocity, determined by the significance of either zonal or meridional or both components with 5% significance, are colored in magenta. The western and eastern magenta boxes represent the Arctic Pacific dipole sea ice index regions.
PIOMAS dataset provides. The formulation of the ocean heat budget can be written as where C p is the specific heat of seawater (3850 J kg 21 K 21 ), r o is the density of seawater (1025 kg m 23 ), T m is the seawater temperature averaged over the control volume, huT*i and hyT*i are the ocean lateral heat fluxes, Q surface is the net heat flux at the surface with positive value downward into the ocean (the Q surface includes both the atmosphere-ocean heat fluxes in the ice-free region and sea ice-to-ocean heat fluxes in sea icecovered region), and R is the residual term including other processes, and x, y, and z denote longitudinal, latitudinal, and vertical directions respectively. In the ocean model producing the PIOMAS product (i.e., the Parallel Ocean Program version 2), when the ocean is at the freezing temperature (i.e., 21.88C) and if additional heat loss occurs due to the air-sea fluxes over the sea ice free region or ocean lateral/vertical fluxes, the ocean temperature stays at 21.88C and the sea ice grows (the interior sea ice temperature can decrease to some extent, but the sea ice bottom temperature is kept at 21.88C). In such a situation, the temperature tendency term [i.e., r o C p ÐÐÐ BS (›T m /›t) dxdydz] on the left-hand side of Eq. (4) equals zero, which is maintained by the balance between the sea ice-to-ocean heat flux (included in the ÐÐ BS Q surface dxdy) and the sum of the other two terms on the right-hand side of Eq. (4). Figure 8a shows the regression of each heat budget term onto the AM BS sea ice index from ON(21) to AM. After JF, the negative temperature tendencies (red line in Fig. 8a) are mainly sustained by the combined cooling effects of lateral heat fluxes (black line in Fig. 8a) and residual term (magenta line in Fig. 8a), while counteracted by the net heat fluxes from the surface (cyan line in Fig. 8a). In AM, the drop of the net surface heat flux anomalies, persistent negative lateral heat fluxes, and enhanced residual together lead to stronger cooling temperature tendency that provides cold ocean conditions, favoring sea ice growth. If we assume the lateral heat flux anomalies below 60 m to be largely consistent with those above the 60-m depth, given the overall barotropic eastward Atlantic inflow and the Norwegian Coastal Current on this continental shelf with maximum bottom depth of ;300 m (Skagseth et al. 2008;Smedsrud et al. 2010), the residual may reflect the upwelling of negative ocean temperature anomalies from below, also a consequence of reduced lateral heat transport. However, we cannot investigate the vertical heat exchange and other processes resulting in the cooling effect of the residual term because 60-m level is the deepest in PIOMAS dataset, and no vertical velocity field is provided.
We next examine the regressed lateral ocean heat flux in the upper 60 m and ocean-to-sea ice heat flux (Fig. 8b). Note that the net surface heat flux into the ocean [Q surface in Fig. 8a and Eq. (4)] is decomposed into the atmosphere-to-ocean heat flux and sea ice-to-ocean heat flux (as explained in section 2a). The latter is then multiplied by 21 to infer ocean-to-sea ice heat FIG. 7. As in Fig. 5, but for the regressions onto AM BS sea ice index. The green box denotes the BS region (408-608E, 688-788N), where the BS sea ice index is defined.
flux anomalies (cyan line in Fig. 8b). The lateral ocean heat flux anomalies are dominated by the zonal heat flux huT*i anomalies, while the meridional component is negligible (not shown). The regressed total zonal heat flux huT*i anomaly through the western face of the BS region at 408E brings less heat into the BS from late autumn to spring, during which the positive BS SIC anomalies gradually increase (black line in Fig. 8b). Furthermore, the evolution of total zonal heat flux anomalies is similar to that of the ocean-to-sea ice heat flux anomalies averaged over the BS (cf. black and cyan lines in Fig. 8b), although the amplitude of latter is smaller. The results suggest that reduced horizonal ocean heat carried into the BS by ocean currents in the upper 60 m alone can account for the reduced amount of heat used to melt the sea ice and cooler ocean conditions that favors sea ice growth above while the atmospheric heat input into the ocean is mostly balanced by the residual, likely the ocean contribution from deeper levels (Fig. 8a). Further decomposition following Eq.
(3) shows that the mean temperature difference carried out by anomalous zonal flow (hu 0 T*i, green line in Fig. 8b) and the anomalous temperature difference carried out by mean zonal flow (huT* 0 i, yellow line in Fig. 8b) each contributes to about 50% of the total zonal heat flux, whereas the anomalous temperature difference carried out by anomalous zonal flow (hu 0 T* 0 i, magenta line in Fig. 8b) contributes little. It is interesting that huT* 0 i and hu 0 T*i play roughly equal roles in producing the BS cooling. This indicates that thermodynamical (temperature change solely) and dynamical (ocean current change solely) processes are equally important. The former could be associated with the relatively colder ocean temperature in the region to the west of the BS, which is produced by atmospheric heat flux loss (see fourth column in Fig. 7). The cooling effect of hu 0 T*i is directly related to the reduced ocean currents, which are likely driven by the atmospheric wind changes, as shown by the regression of JF geopotential height anomalies at 1000 hPa onto the AM BS sea ice index (bottom right panel in Fig. 9). Indeed, the positive geopotential height anomalies extending from southeast of Greenland to the BS favor anomalous easterlies (blowing westward) that weaken the prevailing eastward ocean current entering the BS (possibly due to resultant negative wind stress curl).

d. Dynamical linkage between the two Arctic regions
In this section, we investigate the possible dynamical linkage between the sea ice anomalies in the two Arctic regions in order to understand the underlying mechanism that gives rise to 7-month lagged relationship. Regressions of the sea ice velocity fields in lead-lag conditions onto both sea ice indices do not show any persistent, significant trans-Arctic oceanic pathway (not shown). Therefore, direct sea ice mass transport from the Pacific Arctic sector to the BS via an oceanic pathway seems unlikely. In addition, sea ice reemergence is not found in the BS SIT and SST evolutions (Fig. 7). Instead, the JF regression map of the near-surface geopotential height anomaly (bottom right panel in Fig. 9) sheds light on a possibility that an atmospheric pathway may be involved in linking the sea ice anomalies in the two subArctic regions.
We thus calculate the regressions onto the AM BS sea ice index of the geopotential height (hereafter Z) anomalies at 1000, 500 (representing the middle troposphere), and 100 hPa (representing the lower stratosphere) backward from SO (21) to JF. Figure 9 shows that a positive anomaly appears in the central Arctic at all levels in SO(21) and persists until JF. As winter approaches, the positive Z anomalies at 100 hPa become elongated and are tilted from the East Siberian Sea toward Greenland, while the near-surface Z anomalies become weaker but with their center of action tending to follow the migration of Z anomalies above. In JF, the near-surface Z anomalies in the North Atlantic become stronger as those in the middle troposphere and lower stratosphere, possibly reflecting the downward propagation of stratospheric circulation anomalies in winter (e.g., Baldwin and Dunkerton 1999). The near-surface Z anomaly pattern in JF resembles the negative NAO pattern in the North Atlantic sector. A similar NAO-like tropospheric response to stratospheric circulation anomalies was also identified in previous studies (e.g., Dunkerton 1999, 2001;Thompson et al. 2002;Scaife et al. 2005;Kolstad et al. 2010;Charlton-Perez et al. 2013). The negative NAO-like pattern in JF reduces the prevailing westerlies in the North Atlantic sector, especially over the Nordic seas and BS regions, which leads to reduced ocean heat transport into the BS as discussed in section 3c. Hence, the dynamical coupling between troposphere and stratosphere could provide an atmospheric pathway in linking the two sea ice anomalies.
Previous studies have shown that the stratospheric circulation anomalies can persist from autumn to winter, and could propagate downward to affect tropospheric circulations in late winter (e.g., Thompson and Wallace 1998;Dunkerton 1999, 2001;Sigmond et al. 2013;Kidston et al. 2015;Scaife et al. 2016) and subsequent spring BS sea ice extent (e.g., Smith et al. 2018). To illustrate that such a persistent stratosphere-troposphere signal is involved in linking the two regional sea ice anomalies, we show the daily evolution of the normalized (at each level) polar cap Z anomalies averaged over 08-3608 and 708-908N from lower troposphere (1000 hPa) to middle stratosphere (1 hPa), regressed onto the AM BS sea ice index (Fig. 10a). These anomalies are averaged with a sliding 61-day window to be consistent with bimonthly-averaged fields. Significant positive height anomalies start to emerge in the lower troposphere (from 1000 to ;700 hPa) in AS(21). They then intensify and extend upward into the stratosphere (;10 hPa), slowing down FIG. 9. Geopotential height anomalies at 1000, 500, and 100 hPa from SO(21) to JF regressed onto the AM BS sea ice index. The geopotential height fields are normalized by the spatiotemporal standard deviation at each level before regression. The magenta contour lines indicate that the regression coefficients are significant at the 5% level, while the cyan contour lines inform the field significances. the stratospheric circulation (not shown) and resulting in positive Z anomalies as shown in Fig. 10a. In the stratosphere these anomalies persist well into ON(21) and D(21)J, but the zonally averaged anomalies become weaker in the lower troposphere (from 1000 to ;700 hPa) around ND(21). In D(21)J and JF, the positive polar cap height anomalies tend to propagate downward toward the lower troposphere, and they become weak and insignificant thereafter. The anomalous circulations in zonal-mean metric are consistent with the spatial regression maps in Fig. 9.
The regression of the polar cap Z anomalies onto the SO(21) dipole sea ice index (Fig. 10b) indicates consistent, but generally weaker signals, in particular below 200 hPa from SO(21) to ND(21). This is mainly due to the smaller extent and the larger spatial inhomogeneity of the positive Z anomalies inside the polar cap regressed on the SO(21) dipole sea ice index (not shown) than when they are regressed on the AM BS sea ice index (Fig. 9), which leads to weaker area-averaged values. Nonetheless, the persistent stratospheric circulation anomalies since SO(21) and seemingly downward propagation of the Z anomalies after ND(21) are also found in Fig. 10b, which supports that stratospheric processes are involved in linking the two regional sea ice anomalies.
Many previous studies have suggested that the Arctic sea ice anomalies could enhance the upward propagation of planetary waves into the stratosphere and affect the stratospheric circulation (e.g., Jaiser et al. 2012;Kim et al. 2014;Peings and Magnusdottir 2014;Kretschmer et al. 2016;Nakamura et al. 2016). Others have shown that the planetary waves propagating into the stratosphere can cause a slowdown of the stratospheric polar vortex (e.g., Karpetchko and Nikulin 2004;FIG. 10. (a) Polar cap (averaged over 08-3608 and 708-908N) daily geopotential height anomalies regressed onto the AM BS sea ice index. The 61-day running averages and normalization by the temporal standard deviation at each level are applied before calculating the regression. (b) As in (a), but regressed onto the SO(21) dipole sea ice index. The regression coefficients within the black contour lines are significant at the 5% level, while the cyan contour lines inform the field significance. (c) The SO(21) vertical W z and zonal W x wave activity flux anomalies regressed onto the SO(21) dipole sea ice index. The vertical W z and zonal W x anomalies are averaged over the 688-788N latitudinal band. For better illustration, W z are multiplied by a factor of 100. The black arrows indicate the wave activity flux regressions are significant at the 5% level. The two magenta thick lines denote the Arctic Pacific dipole area. Polvani and Waugh 2004;Kidston et al. 2015). Hence, to investigate if the sea ice dipole anomalies are related to such processes, we examine the SO(21) vertical and horizontal wave activity flux (W x and W z ) anomalies, calculated following Takaya and Nakamura (2001), regressed onto the SO(21) dipole. Significant upward W z anomalies indeed appear above the western part of the Arctic Pacific sea ice dipole region and extend from the surface to the stratosphere, while the W z anomalies present an overall noisy pattern elsewhere (Fig. 10c). Similar results are found in the regression of W z onto the dipole in ON (21), but with weaker strength and less significance as the amplitude of sea ice dipole anomalies becomes smaller in ON(21) (not shown). The near-surface temperature gradient associated with the dipolar sea ice anomalies could be the forcing mechanism, but further dynamical analyses and numerical experiments are needed to confirm the underlying mechanism and establish how these vertical wave fluxes interact with the stratospheric circulation in late autumn and early winter.

Summary and discussion
This study finds that the springtime BS sea ice anomalies are closely related to the autumn Arctic Pacific sea ice dipole anomalies in the previous year. Correlation analysis reveals that ;50% of AM BS sea ice variability can be explained by SO(21) Arctic Pacific sea ice dipole variations, informing a new source of BS sea ice predictability with 7-month lead time. Indeed, a cross-validated linear regression method using the autumn Arctic Pacific sea ice dipole indicates high prediction skills, with 0.54-0.85 ACCs, for the subsequent springtime BS sea ice condition. The result suggests that inclusion of the Arctic Pacific sea ice dipole information at 7-month lead time could improve springtime BS sea ice seasonal forecast skill. The prediction skills are high particularly in subperiods before the mid-1990s and after the mid-2000s, but weaker in subperiods centering on the late 1990s. The reasons why the prediction skills experience nonstationary changes require further investigation and need to be considered in the regional BS sea ice prediction system.
We have performed case studies and regression analysis to understand the physical meaning of autumn sea ice dipole anomalies in the Arctic Pacific sector. Persistent shortwave radiative anomalies from late spring to early autumn, possibly related to cloud cover variability, are shown to be the major driver of the autumn sea ice dipole anomalies, which reflects the studies on the effect of preceding Arctic cloud on autumn sea ice variation (e.g., Kay and Gettelman 2009;Cox et al. 2016;Huang et al. 2019). Meanwhile, the shortwave-albedo feedback could enhance the persistence of the sea ice dipole anomalies. Previous studies suggested that Arctic sea ice variability may be affected by large-scale teleconnections (e.g., Screen and Francis 2016;Zhang et al. 2020). However, we find no evidence that similar teleconnections from outside the Arctic could force the dipole sea ice anomalies, nor do we find persistent concomitant snow and SST anomalies in any other regions of the Northern Hemisphere (even in the tropical Pacific) from MJ(21) to SO(21) (not shown). Therefore, the essential driver of the sea ice dipole, which effectively affects the Arctic cloud cover, does not result from forced largescale teleconnections. On the other hand, the oceanic heat transport from the Bering Strait into the Arctic Ocean may not contribute to the sea ice dipole variability either. We indeed found that the ocean-to-sea ice heat flux anomalies are very small preceding the autumn sea ice dipole. However, the stratospheric circulation anomalies in spring have been shown to be able to affect tropospheric circulation and drive Arctic Pacific sea ice dipole anomalies (Kelleher et al. 2020). As such, whether the sea ice dipole anomalies can be treated as a physical identity requires a further comprehensive understanding on the role of these potential physical processes play.
We further show that atmospheric circulation anomalies extending from the North Atlantic sector to the Nordic seas play important roles in driving lateral ocean heat transport changes that affect the BS sea ice extent. The relationship between BS sea ice extent and lateral ocean heat flux carried out by anomalous ocean currents is consistent with previous findings (e.g., Yeager et al. 2015;Zhang 2015;Årthun et al. 2017;Lien et al. 2017;Oldenburg et al. 2018). Although these studies were aiming for the BS sea ice predictability at longer time scales, it would be interesting to apply a similar approach to understand whether or not inclusion of the ocean dynamics and its variability could improve the BS sea ice condition at seasonal time scale. Our analyses further show that the physical processes linking the two Arctic regions involve an atmospheric pathway that extends the autumn sea ice dipole influence to the winter and leads to a negative NAO-like signal in the troposphere via stratosphere-troposphere coupling. The anomalous tropospheric circulation then reduces the horizontal ocean heat fluxes into the BS region during winter and early spring, resulting in a cooler ocean condition and increased BS sea ice extent. Since the winter NAO seems to play a key role in linking sea ice in the two regions, we also included NAO information in our empirical prediction model (cyan line in Fig. 3a). The prediction skill is then slightly enhanced in most of the 12-yr intervals selected in the analysis period (magenta line in Fig. 3b). It is noted that a previous study also emphasized the roles of winter stratospheretroposphere coupling and subsequent tropospheric NAO anomalies in predicting following spring BS sea ice condition (Smith et al. 2018). Our results support their findings and extend the spring BS sea ice predictability further back to previous autumn.
It has been suggested that the BS sea ice changes may interact with nearby Ural blocking systems to form a positive feedback that reinforces the BS sea ice anomalies (e.g., Luo et al. 2016aLuo et al. ,b, 2018). Based on regressed sea level pressure in early spring (not shown), we indeed find that negative sea level pressure anomalies with a center of action south of Novaya Zemlya Island and north of the Ural Mountains (approximately at 658N, 608E) emerge when the BS sea ice anomalies are intensified after AM. The resultant northerly wind anomalies in the BS region may additionally contribute to the growth of sea ice and the maximum sea ice anomalies found in AM.
Finally, the finding of this study is based on regression and correlation analysis. However, correlation does not necessarily imply causality. Hence the inferred dynamical processes for linking the two SIC anomalies remain somewhat speculative, albeit physically plausible, and they require further analysis. In future studies, we anticipate using advanced statistical tools, such as causal network analysis (Runge et al. 2015), which was applied to explore the causal link between the Arctic Oscillation and Barents-Kara sea ice extent (Kretschmer et al. 2016). A hierarchy of carefully designed regional and global climate model simulations will also be considered to test the dynamical linkage presented in this study.